{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Componenet Analysis (PCA) Multi-Node Multi-GPU (MNMG) Demo\n",
    "\n",
    "PCA is a dimensionality reduction algorithm that works really well for datasets which have correlated columns. It combines the features of X in linear combination such that the new components capture the most information, or variance, of the data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike the single-GPU implementation, The MNMG PCA API currently requires a Dask cuDF Dataframe as input. `transform()` also returns a Dask cuDF Dataframe. The Dask cuDF Dataframe API is very similar to the Dask DataFrame API, but underlying Dataframes are cuDF, rather than Pandas.\n",
    "\n",
    "For information on converting your dataset to Dask cuDF format: https://rapidsai.github.io/projects/cudf/en/stable/dask-cudf.html#multi-gpu-with-dask-cudf\n",
    "\n",
    "For more information about cuML's PCA implementation: https://rapidsai.github.io/projects/cuml/en/stable/api.html#cuml.dask.decomposition.PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "\n",
    "import pandas as pd\n",
    "import cudf as gd\n",
    "\n",
    "from cuml.dask.common import to_dask_df\n",
    "from cuml.dask.datasets import make_blobs\n",
    "\n",
    "from dask.distributed import Client, wait\n",
    "from dask_cuda import LocalCUDACluster\n",
    "\n",
    "from dask_ml.decomposition import PCA as skPCA\n",
    "from cuml.dask.decomposition import PCA as cumlPCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Start Dask Cluster\n",
    "\n",
    "We can use the `LocalCUDACluster` to start a Dask cluster on a single machine with one worker mapped to each GPU. This is called one-process-per-GPU (OPG). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster = LocalCUDACluster(threads_per_worker=1)\n",
    "client = Client(cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 2**18\n",
    "n_features = 20\n",
    "\n",
    "n_components = 2\n",
    "whiten = False\n",
    "random_state = 32\n",
    "svd_solver = \"full\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Data\n",
    "\n",
    "### GPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "X_dcudf, _ = make_blobs(n_samples, \n",
    "                        n_features, \n",
    "                        centers=1, \n",
    "                        cluster_std=0.01,\n",
    "                        random_state=random_state)\n",
    "wait(X_dcudf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Host\n",
    "\n",
    "Dask-ML accepts a Dask.Array, instead of Dask.Dataframe, as input. Dask ML also wants to know the exact sizes of the partitions so we use the argument `lengths=True` to get this information from the workers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ddf = to_dask_df(X_dcudf).to_dask_array(lengths=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scikit-learn Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "pca_sk = skPCA(n_components=n_components,\n",
    "               svd_solver=svd_solver, \n",
    "               whiten=whiten, \n",
    "               random_state=random_state)\n",
    "\n",
    "result_sk = pca_sk.fit_transform(X_ddf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## cuML Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "pca_cuml = cumlPCA(n_components=n_components,\n",
    "                   svd_solver=svd_solver, \n",
    "                   whiten=whiten,\n",
    "                   random_state=random_state)\n",
    "\n",
    "result_cuml = pca_cuml.fit_transform(X_dcudf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate Results\n",
    "\n",
    "### Singular Values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "passed = np.allclose(pca_sk.singular_values_, \n",
    "                     pca_cuml.singular_values_.to_array(), \n",
    "                     atol=0.01)\n",
    "print('compare pca: cuml vs sklearn singular_values_ {}'.format('equal' if passed else 'NOT equal'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explained Variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "passed = np.allclose(pca_sk.explained_variance_, \n",
    "                     pca_cuml.explained_variance_.to_array(), \n",
    "                     atol=1e-6)\n",
    "print('compare pca: cuml vs sklearn explained_variance_ {}'.format('equal' if passed else 'NOT equal'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explained Variance Ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "passed = np.allclose(pca_sk.explained_variance_ratio_, \n",
    "                     pca_cuml.explained_variance_ratio_.to_array(), \n",
    "                     atol=1e-6)\n",
    "print('compare pca: cuml vs sklearn explained_variance_ratio_ {}'.format('equal' if passed else 'NOT equal'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sk_components = np.abs(pca_sk.components_)\n",
    "cuml_components = np.abs(np.asarray(pca_cuml.components_.as_gpu_matrix()))\n",
    "\n",
    "passed = np.allclose(sk_components, cuml_components, atol=1e-3)\n",
    "print('compare pca: cuml vs sklearn components_ {}'.format('equal' if passed else 'NOT equal'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "passed = np.allclose(result_sk, np.asarray(result_cuml.compute().as_gpu_matrix()), atol=1e-1)\n",
    "print('compare pca: cuml vs sklearn transformed results %s'%('equal'if passed else 'NOT equal'))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
